Phenomenological Relations for Axial Quasi-normal Modes of Neutron Stars with 

Realistic Equations of State 
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Here we investigate the axial w quasi-normal modes of neutron stars for several equations of 
state. In particular, we study the influence of the presence of hyperons in the core of the neutron 
stars. We have obtained that w-modes can be used to distinguish between neutron stars with 
hyperons and without hyperons for compact enough stars. We present phenomenological relations 
for the frequency and damping times with the compactness of the neutron star for wl and wll 
modes showing the differences of the stars with hyperons in the core. Also, we obtain a new 
phenomenological relation between the real part and the imaginary part of the frequency of the w 
quasi-normal modes, which can be used to estimate the central pressure of the neutron stars. We are 
able to construct explicitly the low compactness limit configuration for which the fundamental wll 
mode vanishes. Finally, we have studied the influence of changes in the cots- crust transition pressure 
obtaining that it is very small. To obtain these results we have developed a new method based on 
the Exterior Complex Scaling technique with variable angle, appropriate for the treatment of the 
exterior part of the quasi-normal modes, which allow us to impose constringent enough conditions 
to generate pure outgoing quasi-normal modes. For the interior part we use a piece-wise polytrope 
approximation for several equations of state. A complete study of the junction conditions have been 
done. 

PACS numbers: 04.40.Dg, 04.30.-w, 95.30.Sf, 97.60.Jd 



I. INTRODUCTION 

Considerable progress has been made in the last years 
on the gravitational wave detectors development. Large- 
scale interferometric gravitational wave detectors, like 
LIGO, GEO, TAMA and VIRGO, have reached the orig- 
inal design sensitivity, and are currently beginning to op- 
erate with frequencies between 1Hz to 1kHz. The sensi- 
tivity is continuously being enhanced [l[ . First detections 
are expected to happen within the next five years. These 
observations are of huge importance because they can 
be used to perform stringent tests of general relativity. 
But also because gravitational wave detection opens a 
new window to observations of the insights of numerous, 
and probably also new, astrophysical processes. Possible 
sources of gravitational waves are interacting black holes, 
coalescing compact binary systems, stellar collapses and 
pulsars 

Neutron stars are major candidates to detectable grav- 
itational wave sources because of their rich emission spec- 
tra, that lays inside the frequency range of current de- 
tectors. A recent review about gravitational waves from 
neutron stars can be found in Q. The study of the sig- 
nature of the equation of state on gravitational radiation 
from neutron star has been developed in 4-9]. The first 
direct search for the gravitational-wave emission associ- 
ated with oscillations of fundamental quadrupole mode 
excited by a pulsar glitch has been presented in [l(| ■ 

Neutron stars are found inside pulsars, and are thought 
to originate after the collapse of massive star cores. The 



supernova explosion causes violent oscillations of the re- 
sulting compact star. Similarly the coalescence of com- 
pact bodies like white dwarfs leaves behind an oscillating 
neutron star. The resulting excess of energy that causes 
the star to oscillate is radiated in the form of gravita- 
tional radiation [2|]. 

It is well known that although the spectrum of neu- 
tron star oscillations is continuous, general perturbations 
causes the star to ring with concrete oscillation frequen- 
cies that dominates over the rest of possible frequencies 
for a certain period of time after the perturbation. These 
resonant frequencies can be studied introducing the con- 
cept of quasi-normal modes, that is, eigenmodes of os- 
cillation, that, although do not form a complete set of 
functions in which to expand every possible perturba- 
tion evolution, are very useful in the determination of 
the eigenfrequencies for which the star tends to oscillate 
[Ill4l3| . These eigenfrequencies are given by a complex 
number. The real part gives us the oscillation frequency 
of the mode, while the imaginary part gives us the inverse 
of the damping time. The quasi-normal mode spectrum 
is quite dependent on the properties of the star, i.e. the 
equation of state. 

Neutron stars are compact objects of very high den- 
sity. Inside of them matter is found at extreme densi- 
ties (10 15 g j cm 3 ) . Current theories predict a layer struc- 
ture for the neutron star, essentially composed by the 
core and the crust. The properties of the crust are very 
different from those of the core, and it is thought that 
this region has a solid crystalline structure similar to a 
metal. Although the high density matter inside the star 
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has exotic properties, the resulting relativistic fluid that 
compose the neutron star can be described as a perfect 
fluid. What we need to know is the equation of state for 
the neutron star matter, but beyond the nuclear densi- 
ties this relation is not well understood. The inner part 
of the nucleus is very model dependent, essentially be- 
cause different populations of particle states may appear 
at those densities. The composition is not exactly known 
and several hypothesis exist [H], UH . Along the core of 
the star and specially in the core- crust interface, first 
order phase transitions are expected to be found in re- 
alistic equations of state. These phase transitions result 
in small discontinuities in the energy density of the star 
matter [ilfl6l|. 

In order to systematize the study of constraints placed 
by astrophysical observations on the nature of neutron- 
star matter, several parameterizations of high-density 
equation of state (EOS) have been introduced: a piece- 
wise polytropic approximation by Jocelyn et al. |17| and 
Lindblom spectral decomposition [l8j]. The advantage 
of these parameterizations is that they reduce the total 
amounts of parameters which modelize the equation of 
state to a more tractable quantity, allowing to constrain 
with observational data the values of the parameters, and 
also usually they enhance a numerical integration scheme 
allowing better precision. Let us note that a very large 
family of equations of state can be fitted to these param- 
eterizations with a very good precision. 

The detection of gravitational radiation from neutron 
stars can be a very useful tool to determine the neutron 
star equation of state. The extraction of the spectrum of 
quasi-normal modes (frequency and damping time) from 
these signals could give important constrains to the equa- 
tion of state, and in consequence, information about the 
behavior of matter at densities beyond nuclear matter 

PH. 

The theoretical study of the quasi-normal mode spec- 
trum considering different models of neutron star with 
different equations of state is then well justified. The nec- 
essary formalism were developed first for quasi-normal 
modes of black holes by Regge and Wheeler [26| and 
by Zerilli [27] . Quasi- normal modes can be differenti- 
ated into polar and axial modes. In these papers it 
is found that the equation describing the quasi-normal 
mode perturbation of the Schwarzschild metric is essen- 
tially a Schrodinger like equation: the Regge- Wheeler 
equation for axial perturbations and the Zerilli equation 
for polar ones. For black holes, both types of modes are 
space-time modes. The formalism was studied in the con- 
text of neutron stars first by Thorne [28l432| . Lindblom 
[33l HH, and then reformulated by Chandrasekhar and 
Ferrari [35U371 ] . and Kojima |38j| . In neutron stars, ax- 
ial modes are purely space-time modes of oscillation (w 
modes), while polar modes can be coupled to fluid oscil- 
lations (although a branch of w modes can be found also 
in polar oscillations) . In this paper we will consider only 
axial modes of oscillation. 

The study of the quasi-normal modes spectrum is com- 



plicated because of several reasons. Quasi-normal modes 
can only be studied numerically. No analytical solu- 
tions are known for physically acceptable configurations 
of neutron stars. Another reason is that quasi-normal 
modes arc found as isolated points scattered on the com- 
plex plane, so that usually an exhaustive scan of this 
plane is necessary in order to find the complete spec- 
trum. Also the numerical study is hindered by the very 
definition of quasi-normal mode itself, that as we will see 
explicitly in the following sections, give rise to diverging 
functions that oscillate an infinite number of times to- 
wards spatial infinity. These functions are not well han- 
dled numerically. 

Several methods have been developed to deal with 
these and other difficulties. For a complete review on 
the methods see the review by Kokkotas and Schmidt 
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Chandrasekhar and Ferrari [35H37J used a slowly 
damped approximation <C which reduces 

the system of equations needed to resolve the problem, 
so that the damping time of the quasi-normal mode can 
be given in terms of the real part. This method was also 
used in the determination of w-modes on polar pertur- 
bations of constant density neutron stars by Kojima et 
al 39]. It is only valid for slowly damped modes, which 
are present in highly compact stars, and is not appropri- 
ate for realistic configurations, although it can be used 
to obtain a first approximation of the fundamental fre- 
quencies. 

In an interesting paper [Zo[ , Kokkotas studied the axial 
spectrum integrating the exterior solution up to a finite 
radius, which allows to impose the outgoing wave behav- 
ior to the exterior solution. The matching between the 
interior and exterior solutions results in the construction 
of a Wronskian at the matching surface, which must be 
zero when both solutions correspond to a quasi-normal 
mode. 

The WKB method [4l[ can be used to approximate the 
outer solutions provided that no back-scattering is found 
in the exterior region. This condition usually is satisfied 
when the imaginary part of the mode is small compared 
with the real part [42| . 

The phase function can be studied in order to deal 
with the oscillatory nature of the perturbation function. 
Studied for black holes by Chandrasekhar and Detweiller 
[43| . the phase function plays a fundamental role in the 
numerical approach made by Andersson et al in 42| and 
in [4J| for neutron stars. The phase is usually a well 
behaved function if there is not back-scattering contam- 
ination in the numerical solution. 

The divergence of the outgoing wave solution of the 
perturbation equations makes the avoidance of incom- 
ing wave contamination a difficult task. This problem 
was dealt for black holes by Andersson in [45| by rotat- 
ing the radial variable into a complex variable parallel to 
the anti-Stokes line. This ap pro ach were used for con- 
stant density neutron stars in 42 1 and [44j together with 
the phase function approach. More recently, Samuelsson 
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et al [46(| used a similar complex-radius approach for a 
constant density configuration which dealt with the di- 
vergence at infinity by integrating the Bondi- Sachs phase 
along a fixed path for the complexified radius parallel to 
the anti-Stokes line. 

Another approach is based on the Leaver continued 
fraction method [43], which allows to impose the outgo- 
ing wave behavior on the border of the star as a self- 
consistent equation that can be satisfied iteratively [48| . 
This method has been used successfully with realistic 
models of neutron stars [l9|, [25| 

As commented above, for fluid modes there is also 
a branch of w-modes. As these are purely space-time 
modes, they do not couple to matter, and the inverse 
Cowling approximation can be used to obtain this part 
of the spectrum [44[ . 

In this work we present a new approach to calculate 
quasi-normal modes of realistic neutron stars. We make 
use of several well-known techniques, like the use of the 
phase for the exterior solution and the use of a complex- 
ified coordinate to deal with the divergence of the out- 
going wave. We also introduce some new techniques not 
used before in this context: freedom in the angle of the 
exterior complex path of integration, a determinant al- 
gorithm to calculate the mode based on general junction 
conditions on the boundary of the star, implementation 
of 34 realistic equations of state in a piece-wise poly- 
trope approximation, use of Colsys to integrate all the 
system of equations at once with proper boundary and 
junction conditions and possibility of implementation of 
phase transition discontinuities. These new techniques 
allow us to enhance precision, to obtain more modes in 
shorter times, and also to study several realistic equations 
of state, comparing results for different compositions. 

In section II we will start presenting the well-known 
general formalism of quasi-normal modes. In section III 
we will present an study of the junction conditions for 
a general matching that can include first order phase 
transitions or surface energy layers. In most of the re- 
sults of the paper we will make use of the junction con- 
ditions without phase transitions or energy layers. In 
section IV we describe the numerical method we have 
used to obtain our results, testing the algorithm to ob- 
tain already known modes for a simple model . In sec- 
tion V we present our results, focusing on four differ- 
ent equations of state: High density SLy (pure nuclear 
matter), ALF4 (mixed quark- nuclear matter), GNH3 and 
BGNIH1 (mixed hyperon- nuclear matter). In section VI 
we finish the paper with a summary of the main points 
of the work. 



II. QUASI-NORMAL MODE FORMALISM 

In order to fix the notation we will present a brief 
review of the formalism used to describe axial quasi- 
normal modes. In the following we use geometrized units 
(c=l,G=l). 



We consider an static spherical space-time with metric 
ds 2 = e 2v dt 2 -e 2X dr 2 -r 2 (d9 2 +sin 2 6dip 2 ). Following the 
original papers (28[- HH, we make perturbations over this 
metric, expanding the perturbation in tensor harmonics, 
and taking into account only the axial perturbations, i.e. 
those which under a parity transformation 8 — > it — 8, 
ip —> 7r + (/>, transform like (— 1)' +1 . Once the free gauge 
transformation is selected (we choose the Reege- Wheeler 
gauge \XM)i the components of the metric perturbation 
are: 



h t e = -h! 



d t Y lri 



, Web = h? m smOdthYir 



i»™"sin( 

We = -h] m ^j-,h r4> = h] m sinfld^™, (1) 

where the angular dependence is explicitly given in terms 
of the spherical harmonics Yi m . In general hf m and h] m 
are functions of t and r. 

The matter inside of the star is considered to be a per- 
fect fluid with stress-energy tensor T^ v = (p + p) u fJ, u L ' — 
pg^ v , where p is the pressure, p is the energy density and 
u is the 4-velocity. Axial oscillations do not modify the 
energy density or the pressure of the fluid. Hence, for ax- 
ial modes, inside of the star only the perturbations on the 
4-velocity of the matter have to be taken into account. 
These perturbations are calculated considering the La- 
grangian displacement vector £* of the fluid, and each 
component of the spatial vector can also be expanded in 
terms of tensor harmonics. So, the Lagrangian displace- 
ment that can couple to the axial perturbations can be 
written as: 



(p+p)r 



sin u 

d e Ytr. 



(2) 
(3) 



where Wi m are perturbative functions of t and r, associ- 
ated with the axial Lagrangian displacement. The Ein- 
stein equation G^ v = 87rT M „ are solved considering the 
terms of same perturbative order and with same angular 
dependence. 

At zero order the equations obtained are well known. 
Inside of the star we have: 



dm 2 
—— = Ai:r p, 
dr 



dp 
dr 



<P+P) 



4irr 3 p 



(r — 1m)r 
dv 1 dp 

dr p + p dr' 



(4) 
(5) 
(6) 



and outside we have the Schwarzschild solution with 
gravitational mass M. 

The equations for first order perturbation are obtained 
in a similar way, but they are time dependent and depend 
linearly in the perturbative functions. Hence, we can ex- 
tract the time dependence from the functions by making 
a Laplace transformation: 



(7) 
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This transformation eliminates the time derivatives but 
introduces explicitly into the equations the eigen- value u>. 
After some algebraic work with the resulting equations, 
it can be seen that the axial perturbations must satisfy 
the well-known Regge- Wheeler equation 27]. 



'2 ylm 



d z z 



+ [to 2 - V(r)] Z lm = 



where r* is the tortoise coordinate: 



dr. 



(8) 



(9) 



and the eigen-frequency of the axial mode is a complex 
number u> = cjsr + iu>Q. The equation is valid both inside 
the star and outside of it, and the potential can be written 
as: 



.,2// 



V(r) 



+ l)r + Airr 3 (p + p) - 6m] 



(10) 



where I is the spherical harmonic index. In this work 
we will consider only the I = 2 case. Inside the star the 
pressure p, and the mass m are functions of r obtained by 
solving equations (J4j) -(|6]), and an equation of state must 
be provided. Outside of the star, p = p = and m = M, 
the gravitational mass of the Schwarzschild metric. 

Note that in general the perturbation function Z is a 
complex function. These translates into a set of two real 
second order differential equations (one for the real part 
Zsft, and another for the imaginary part Zq). 

Z has to satisfy a set of boundary conditions that can 
be obtained from the following two requirements [351 ] : i) 
the perturbation must be regular at the center of the star, 
and ii) the resulting quasi-normal mode must be a pure 
outgoing wave. 

The first requirement imposes the following Taylor ex- 
pansion to the Z function at the origin 



Z = Z^r l+1 +o(r l + 3 ) 



(11) 



Note that Z^ is a non-fixed parameter, and we will make 
use of this freedom in our numerical method of quasi- 
normal mode calculation. 

The second requirement imposes the quasi-normal 
modes to behave as purely outgoing waves at radial infin- 
ity. In general a quasi-normal mode will be a composition 
of incoming and outgoing waves, i.e. 



lim Z out ~ e~ lujr * . 



(12) 
(13) 



Note that, while the real part of ui determines the oscil- 
lation frequency of the wave, the imaginary part of the 
eigen-value determines the asymptotic behavior of the 
quasi-normal mode: If we call r = l/u>i, the ingoing 
mode will behave as: 



lim Z T 



,-r„/T 



(14) 



which is a convergent function, and the outgoing mode: 

(15) 



lim Z out ~e r + /T , 



which is a divergent function. Outgoing quasi-normal 
modes are divergent at radial infinity, while ingoing ones 
tends exponentially to zero as the radius grows. Hence, 
the purely outgoing condition does not give us a cor- 
rect boundary condition that could be satisfied at radial 
infinity. In principle, the purely outgoing quasi-normal 
mode condition could be imposed at a very distant point, 
but because of the quick convergence of the ingoing sig- 
nal, which tends exponentially to zero at radial infinity, 
every small numerical error in the imposition of this be- 
havior will be amplified as we approach the border of 
the star, resulting in a mixture of outgoing with ingoing 
waves. The completely different asymptotic behavior of 
the outgoing and ingoing quasi-normal modes makes very 
difficult to obtain quasi-normal modes without any mix- 
ture of ingoing waves. So trying to impose this behavior 
numerically is probably not the best treatment because 
small numerical deviations far away from the star will 
translate into big incoming quasi-normal modes contri- 
butions at the surface of the star. This is precisely what 
happens when attempting to calculate strongly damped 
modes. 

We will deal with this problem as we will see below. 



III. JUNCTION CONDITIONS 

To determine completely the problem, we must study 
the junction conditions between the interior solution and 
the exterior solution on the boundary of the star. The 
junction conditions between two space-times have been 
considered in classical works by Darmois [49| , Lichnerow- 
icz [SO], O'Brien and Synge [5l|, and Israel [52}. Junction 
conditions for general hypersurfaces have been considered 
in [53J , Symmetry-preserving matchings have been stud- 
ied in |54j ]. and an analysis of how the induced metric 
and extrinsic curvature of a hypersurface change, when 
space-time metric is perturbed to second order and the 
hypersurface itself is deformed perturbatively to second 
order, is presented in [55j. Junction conditions for rel- 
ativistic rotating stars have been considered in [56| and 
57} . The junction conditions in surface of energy discon- 



tinuity have been considered in [58] to study the influence 
of phase transition in neutron star quasi-normal modes. 

We will impose the usual junction conditions. But, 
in order to have the possibility to study the influence of 
small changes in the core- crust transition pressure, we 
will include the more general case in which the star is 
surrounded by a surface layer of energy density. 

We use Darmois conditions, the intrinsic formulation 
of these junction conditions, which impose certain con- 
strains in the continuity of the fundamental forms of the 
matching hypersurface S. These conditions are the fol- 
lowing: the first fundamental form has to be continuous 
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across S, and the jump in the second fundamental form 
will be given in terms of the stress-energy tensor of the 
surface layer 

This surface S is defined by the points where the pres- 
sure is null, or more generally, constant. As we have 
seen previously the axial oscillations do not modify the 
pressure or the density, so the surfaces of constant pres- 
sure of the perturbed star are essentially the same as the 
surfaces of constant pressure of the static star: 



p = p^ (r) = const. 



(16) 



This condition is equivalent to requiring r — const. = R. 
We calculate a vector normal to these surfaces, n^, and 
from it the projector h uv = + n^n^ and the pro- 
jected covariant derivative of the normal vector = 
h"hP,n a; p. Evaluating these tensors at r = const. = R, 
from inside and from outside S, we obtain the first fun- 
damental form and the second fundamental form inside 
and outside S [5!| . We must make a series expansion of 
these conditions in order to obtain matching relations for 
each one of the metric functions. We also must expand in 
spherical harmonics and make a Laplace transformation, 
in order to obtain the junction of the radial component 
of the functions. 

Let us consider a surface S of constant pressure inside 
the star, and a surface layer on top of S described by 
a surface stress-energy tensor of a perfect fluid of the 
following form: 



Tg»(R)=su»(R)u» s (R), 



(17) 



where, S is defined by the points where r — const. = R, 
e is the surface energy density of the fluid moving in 
S, and Ug(R) is its velocity. The continuity of the first 
fundamental form gives us two conditions: the continuity 
of v 



Vext, 



and the continuity of perturbation function Z: 



ylm rylra 

^int ~ ^exf 



(18) 



(19) 



The second fundamental form conditions impose a 
jump in the pressure if e ^ 0. Let the pressure in the 
inner part of the surface S, be Pinti an d in the outer part 
be 

Pext 

with p ext < Pint- 
Making the same steps as with the first fundamental 
form, we obtain the following conditions for the zero or- 
der functions: 



M ext = M int + 1 



2Af„ 



R 



e - %n z R A e\ (20) 



M ext + 4irR 3 p ext M mt + 4ttR 3 Pv 



R\ll-™^ 



R 2J 1 _2M^ 



Aire. (21) 



In this work we will consider p ext — 0. That is, the 
surface of star with pressure p = pint is being matched 



directly to the empty exterior surface of the star, in which 
case, M ext is just the Schwarzschild mass M. If we choose 
Pint as the core- crust transition pressure, this matching 
condition is equivalent to approximating exterior layers 
of the star to a surface layer energy density that envelops 
the core of the star (surface crust approximation). The 
equation (|20l) can be interpreted as giving the total mass 
of the star to zero order (M ext ), as the addition of three 
terms, the first one representing the core mass (Mj„ t = 
4-7T J Q pr 2 dr), the second one representing the mass of 
the crust and the third one being a negative bounding 
energy term. The radius of the resulting configuration 
will be the radius of the core of the star. Equation (|2"Tj) 
can be seen as a condition for the determination of the 
radius of the core of the star when approximating the 
crust of the star to a thin surface energy density. 

In the surface crust approximation, these relations al- 
low us to use the core-crust transition pressure Pi n t, as a 
new parameter of the resulting configuration. 

Up to first order in the perturbative expansion we ob- 
tain junction conditions for the axial functions. Essen- 
tially we obtain the following condition for the derivative 
of the Z function: 



dZ l 



dZ L 



dr 



dr 



R 



(22) 



Together with condition (|19p we have enough information 
to perform correctly the matching of the perturbation. 

Let us note that, the usual treatment without surface 
crust approximation is recovered by making in the previ- 
ous conditions p^t = Pext = and e = 0, in which case 
the continuity of the metric functions and of their deriva- 
tives is obtained. This is what we do in almost all the 
paper except in the final part of the results section, where 
the influence of the changes in the core-crust transition 
pressure is analyzed. 

Because the junction conditions for the axial pertur- 
bations (jT9"]l and (|2"21 are homogeneous in Z, we have a 
freedom in the global multiplicative factor (the free co- 
efficient from equation (TTTT) ). Back to the problem of 
the determination of the purely outgoing quasi-normal 
modes, note that now we have sufficient conditions in or- 
der to solve the problem: regularity at the origin for the 
interior perturbation, asymptotic behavior for the exte- 
rior one, and matching conditions for both the function 
and the first derivative of the perturbation. 



IV. NUMERICAL METHOD 

In this section we will explain the numerical method 
we use to obtain outgoing quasi-normal modes of realistic 
stars. We implement this method into different Fortran 
programs and routines, making use of the Colsys pack- 
age [60l ] to solve numerically the differential equations. 
The advantage of this package is that it allows the uti- 
lization of quite flexible multi-boundary conditions, that 
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as we will see in the following sections, are important 
in our method of quasi-normal mode determination. On 
the other hand, Colsys allows to construct an adaptative 
mesh of points depending of the precision we would like 
to achieve. In this manner, the more precision we need on 
the perturbation function, the more points Colsys intro- 
duces into the mesh. As we will see, the static solution 
is calculated at the same time as the perturbation, so 
the static solution functions will be also adapted to the 
precision we require to the Regge- Wheeler function. Col- 
sys also allows the use of the continuation method: we 
can use as an initial guess for a new numerical solution 
another numerical solution which parameters are close 
enough to the parameters of the new solution. 

We will start explaining how the outgoing modes are 
determined in the exterior region of the metric. 



A. Calculation of outgoing quasi-normal modes 
outside the star 

As it has been introduced in previous sections, quasi- 
normal modes has two main numerical problem. Quasi- 
normal modes present an oscillating part given by the 
real part of the eigen-frequency, that remains towards 
infinity, as can be seen in relation (j 13[) . Also, quasi- 
normal modes diverge towards spatial infinity because of 
their pure outgoing wave behavior. 

In order to deal with these two main issues we will 
follow HH-EH, using the phase function together with a 
rotation of the radial coordinate into the complex plane. 

The phase function (logarithmic derivative of the Z 
function) , 



g = Y > 



(23) 



does not oscillate towards asymptotic infinity. Even 
more, note that the asymptotic behavior of the phase, 
considering a general combination of outgoing and ingo- 
ing waves is: 



g — > iuj 



Ae 1 ' 



Be- 



rn 



When r is real, there are several possibilities. If the mode 
is purely ingoing, B=0, then g — > iuj. If the mode is 
purely outgoing, A = 0, or if it is a mixture, A ^ 0, then 
g — > —iuj. Although the phase allows us to eliminate 
the oscillation, we can not distinguish between the mixed 
incoming and outgoing waves and the purely outgoing 
ones. The equation for the phase can be written as: 

d 9 2 , u 2 r 2 , 2M - l(l+l)r-6M _ p r9 ,-s 
+ ( 2 -2M)i + r(r-2M)9 ~ r i(r-2M) ~ U ' ( Zi> ) 

As we will see, the phase function has all the informa- 
tion we need to obtain the quasi-normal mode, provided 
enough junction conditions. 



The divergence problem can be avoided by rotating the 
radial coordinate into the complex plane. This technique 
of complexification of the integration variable is called 
Exterior Complex Scaling |6l| - t63| . An analytical exten- 
sion of the g function, which is a function of r, is made 
into a g function that is dependent of the generalized 
complex coordinate. This new function has no physical 
meaning, but the eigen-values obtained by integrating 
the equation along the extended radius will remain unal- 
tered. 

We parameterize the complex coordinate along a line 
in the complex plane. We will take a straight line with 
an angle a with respect the real axis. 



R + ye~ 



(26) 



with y G [0,oo). The g function has the same limits in 
the case of a purely outgoing or a purely ingoing wave 
than with the real coordinate. But the mixed waves now 
depend on the angle of the path of integration a. It 
can be demonstrated that a mixed ingoing-outgoing wave 
now presents the following asymptotic behavior for its 
phase: 

wsr sin a + ojq cos a > g — > —iuj, (27) 
wsr sin a + ojq cos a < =>- g — >■ iuj. (28) 

If wso sin a + ujn. cos a — 0, that is, when a = — arctan — 
then the phase become an oscillatory function. This limit 
line L is parallel to the anti-Stokes line. This path is not 
useful as our objective is to compactify in order to con- 
trol completely the boundary conditions. The interesting 
part of the plane is beyond the line L, because in this re- 
gion the asymptotic behavior of the purely ingoing waves 
are exactly the same as the asymptotic behavior of the 
mixed waves, g — > iuj. In other words, beyond the line L 
the asymptotic behavior of the outgoing waves g — > —iuj, 
prohibits any ingoing wave contribution. So now, if we 
compactify the parameter y by defining x = r^j, we 
can impose at infinity (x — 1), the boundary condition 
g\ x =i = —iw. This boundary condition impose to the 
obtained solution a purely outgoing wave behavior. 

In our method, the angle of the integration path is 
treated as a free parameter of the solution. This an- 
gle can be adjusted appropriately to increase precision 
and integration time, but always remaining in the region 
where g — > —iu. Exterior Complex Scaling with variable 
angle is widely used in other contexts like atomic and 
molecular physics, where it gives excellent results [64|. 

The use of the phase and Exterior Complex Scaling 
with variable angle allows us to compactify, so we can 
impose the outgoing quasi-normal mode behavior as a 
boundary condition at infinity without using any cutoff 
for the radial coordinate. 

Note that the equation we must solve now is a first 
order equation for a complex function. So we have two 
first order differential equations for two real functions, 
gsjf and <7q, together with these two boundary conditions: 
g${\ x= i = ojq and gQ\ x =i — — ^sr- So we can obtain 
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a phase function given a pair of values (wsr, wg). We 
still need the interior part of the quasi-normal mode to 
determine if this pair of values corresponds to an eigen- 
frequency or not. 



B. Calculation of quasi- normal modes inside a 
realistic star 

Inside the star the Regge- Wheeler equation (JS) can be 
solved without any problem. The perturbation function 
is oscillating, but the number of oscillations will be Unite 
between r = and r = R. 

To obtain the perturbation function, we need to solve 
the equations for the zero order functions (j3|, © and (JB]). 
To do so, we integrate them with the proper boundary 
conditions: regularity at the origin and junction condi- 
tions (PHI , using a realistic equation of state. 

In order to have the border of the star perfectly under 
control, we introduce the radius of the star as a function 
with its own differential equation (constant everywhere 
inside the star) into the complete system of equations. 
In principle the interior of the star is a free boundary 
problem, where the border of the star is obtained where 
the pressure vanishes. Using the radius as an unknown 
function allow us to integrate the configuration imposing 
physical conditions on the boundaries, and obtain the ra- 
dius as a result in the same manner as one obtains the 
mass. Colsys is the appropriated way of solving the prob- 
lem, because Colsys allows to solve differential equations 
with quite general boundary conditions. 

To implement the equation of state we use a piece- wise 
polytrope approximation, done by Read et al [17| . In this 
approximation the equation of state is approximated by a 
polytrope in different density-pressure intervals. For den- 
sities below the nuclear density, the SLy equation of state 
is considered. For higher densities three different regions 
are used. Then six parameters are needed to model the 
equation of state: two dividing densities, three adiabatic 
indices and a value of the pressure at an endpoint of one 
of the intervals. The piece- wise polytrope approximation 
is done in [l7j for 34 different equations of state, mini- 
mizing the error. We have implemented into our routines 
this approximation for all these equations of state so that 
once it is selected, the programs read the corresponding 
6 parameters and Colsys begins the integration, applying 
each polytrope to the corresponding density intervals. 

As we commented previously, Colsys uses an adapta- 
tive mesh of points to generate the solutions, so that if 
more precision is needed in a certain region of integration, 
Colsys automatically generates more points there. Be- 
cause the Regge- Wheeler equation depends on the pres- 
sure and mass, this feature is very interesting because al- 
lows the Z function to have more points in specially com- 
plicated regions where, for example, the pressure drops 
quickly to zero, increasing the precision in the calcula- 
tion of the quasi-normal mode. To take advantage of 
this, we decide to integrate together the zero order and 



the Regge- Wheeler equation. In this way, we allow Col- 
sys to generate more points in the zero order functions 
if we require a more precise result in our quasi-normal 
modes. Colsys allows the use of different tolerances for 
each function, so different precisions can be achieve in the 
zero order or in the perturbation functions if necessary. 

The Regge- Wheeler equation is a second order equa- 
tion for a complex function Z. That is, we need to solve 
two second order equations for two real functions Zsr and 
Zcj. We need two boundary conditions for each one of 
the functions. 

We must impose the regularity condition at the origin. 
Numerically the condition Z = is not very constringent, 
so we choose to redefine the function. Because we are in 
1=2 modes, we define 

Z = Z/r 3 , (29) 

so that we integrate the equation for Z. We choose the 
following two boundary conditions at the origin 

Z' = 0. (30) 

These conditions impose a regular behavior at the ori- 
gin. We still need another two conditions. The remaining 
boundary conditions must guarantee somehow that the 
matching between the interior and the exterior solutions 
correspond to an outgoing quasi-normal mode. To deter- 
mine such solution, we use what we call the determinant 
method. 



C. Determination of quasi- normal modes using the 
Junction Conditions 

In order to test if for a certain couple of values (wsr , ug) 
the obtained solutions are a quasi-normal mode, we must 
check that the boundary conditions (fH?|) and (|2"2")l are 
satisfied. We use these junction conditions to determine 
the quasi-normal using the following method: 

A general solution to equation ©, once the 2 regularity 
conditions are imposed, can be obtained by generating to 
independent solutions and combining them. Note that 
we still have the free scale factor, so we can obtain the 
independent solutions by generating two solutions with 
different internal phases i.e. different Zsr/Zq. relation. 
For example we generate a solution A with Z A \ r= R = 1 
and Z^\ r= R = 2 and another solution B with Z B \ r= R = 
2 and Z§\ r =R = 1. The obtained solutions Z A and Z B 
could be combined to generate any solution Z = CaZ a + 
CbZ b , where Ca and Cb are two constants. 

From the integration of each particular solution we will 
obtain the value of the derivative at the border of the 
star, ^r\ r =R and ^r\ r =R- These values of the Regge- 
Wheeler function correspond to the interior values of the 
matching conditions (those obtained in the interior face 
of the matching surface at r = R). Let us rewrite these 
in the following form, to make it explicit where they are 
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evaluated, if in the outer face of the surface (+) or the 
inner part (-): 



dr 



-R = Z A {-) 
dZ A 



\r=R 



dr 



R, 



(31) 
(32) 



and the same for the solution B. 

From the boundary condition (|19p we know that on 
the outside face the Regge- Wheeler function must have 
the same values as the solution on the inside. Outside 
the star we have calculated the phase, so no value on the 
surface of the star is imposed to the exterior part of the 
quasi-normal mode. We use the junction condition (fT9|) 
to obtain that value. So we have: 



Z A (+) = Z A (-) = 1 + 2*, 
Z B {+) = Z B (-)=2 + i. 



(33) 
(34) 



The values of the derivatives of the Regge- Wheeler func- 
tions on the inner face of the surface are obtained numeri- 
cally from the integration of the Regge- Wheeler equation 
inside the star: ^f— (— ) and ^r( — )• But we also need 
the value of the derivative of the Regge- Wheeler function 
outside the star. This is easy to obtain from the phase 
function g outside the star. From the integration of the 
phase of the Regge- Wheeler function outside the star we 
have g\ x =o — <?(+), the phase evaluated on the outer face 
of the matching surface. By doing the following calcula- 
tion we obtain the exterior derivative: 

r\7 A A7 B 

' Z -(+) = Z A (+)g(+), 5£_(+) = Z B (+)g(+). (35) 



dr 



dr 



This step can be done because we are imposing the con- 
tinuity of the Regge- Wheeler function for each one of the 
solutions calculated inside the star, solutions A and B. 
But this can be perfectly achieved because we still have 
the freedom in the scale factor for the Regge- Wheeler 
function outside of the star. 

Now that the exterior derivative is calculated, we have 
all the necessary data to evaluate the junction condition 
(|22l) (The necessary values of the zero order functions can 
be evaluated using the zero order junction conditions). 
Continuity of particular solutions A and B grants the 
continuity of the composed solution. The matching (|2"2")l 
must be done on the composed solution. The condition 
can be written as: 



dr 



dr 



a A(H 



=A(-) 



dr 



dr 



1 

R 



[C A Z A (+) + C B Z B (+)] =0.(36) 



This condition can be rewritten in matrix form. Note 
that the Regge- Wheeler function is a complex number. 



Both the real part and the imaginary part of relation (|36|) 
must be null. In other words: 



M 



C A 
C b 



(37) 



with M a matrix which components could be written as: 



M] = $(+)- 



,M+) dZ) 



H- 



Z H-) 

R 



,xi- 



1 



(38) 



where i = {A, B} and j = {5ft, 5}. Because constants 
Ca and Cb can not be zero at the same time, the quasi- 
normal mode condition becomes 



det(M) = 0. 



(39) 



The determinant can be calculated numerically from the 
integrated solutions. For each couple of values (wsr,wcj) 
we can generate interior and exterior solutions, and then 
calculate the determinant. In principle it will be dis- 
tinct from zero, but when wsr and wq correspond to a 
quasi-normal mode eigen-frequency, the determinant will 
be zero. Because the boundary conditions and junction 
conditions are only accomplished when we have a outgo- 
ing quasi-normal mode, the determinant will be zero if 
and only if the eigen- value u> corresponds to an outgoing 
quasi-normal mode. 

An important new feature is the use of the continua- 
tion method, which allows to use as initial guess solution 
a previously generated numerical solution. This trans- 
lates into better convergence times and precision. We 
start guessing with some values for o>sr and ujq, and then 
we start moving the parameters, using the numerical so- 
lution as an initial guess to generate the next solution 
with the new parameters. This procedure allows us to 
sweep the complex plane in reasonable times. 

So, in a summary, the method is the following one: 

First we select an equation of state from the 34 equa- 
tions implemented from (l7| . Choosing a central pres- 
sure, a transition pressure and a value of wsr and wg, 
we integrate the static functions with the appropriate 
boundary conditions together with the Regge- Wheeler 
function inside of the star, to which we impose regularity 
at the origin and two different set of boundary conditions 
at the border of the star. We obtain two independent so- 
lutions for Z. Then we integrate the exterior solution (the 
phase function) over a rotated spatial coordinate (which 
is compactified) with an angle beyond the anti-stokes line 
(this angle can be changed in order to optimize the con- 
vergence time and precision of the phase). Imposing the 
outgoing quasi-normal mode condition, which in this co- 
ordinate is well defined, we obtain an outgoing phase 
function. Now, reading the numerical values from the 
results, we evaluate the determinant. 

We must explore the plane (wsRjO-'s) looking for null 
determinants. This process can be automated so that 
the program looks for minimums of the determinant on 
the plane. Colsys needs an initial guess solution in order 
to start the integration algorithm. To generate the first 



9 



TABLE I: Frequency and damping time of the trapped modes 
(t), the spatial modes and the interface modes (i), in central 
density units, for axial 1=2 modes, 2M/R = 0.885. 



^3t 
















(t)0.21386387 2.4320459 


io- 


9 


1.68347876 


6.0243716 


10" 


2 


0)0.29101153 7.7472523 


10" 


-8 


1.76704836 


6.1930099 


10" 


2 


(t)0.36799864 1.0725279 


10" 


6 


1.85066583 


6.3390024 


10" 


2 


(t)0.44463493 9.5187522 


10" 


6 


1.93429817 


6.4699961 


10" 


2 






10" 


5 


2.01793447 


6.5897686 


10" 


2 


0.59557900 


3.3806006 


io- 


4 


2.10157798 


6.7019361 


10" 


2 


0.66896067 


1.4322666 


io- 


3 


2.18522398 


6.8080510 


io- 


2 


0.74090953 


4.5080077 


io- 


3 


2.26888071 


6.9098878 


io- 


2 


0.81280272 


1.0135742 


io- 


2 


2.35254774 


7.0076042 


io- 


2 


0.88590397 


1.7259264 


10" 


2 


2.43622599 


7.1018102 


10" 


2 


0.96037307 


2.4480195 


10" 


2 


(i)l. 45115909 


5.4011419 


10" 


2 


1.03608110 


3.1027484 


io- 


2 


(i) 1.73433905 


1.3945066 


10" 


2 


1.11308152 


3.6685221 


10" 


2 


(i)2. 01106806 


2.2902687 


10" 


2 


1.19147020 


4.1559313 


10" 


2 


(i)2.29393679 


3.1975625 


10" 


2 


1.27122330 


4.5838031 


io- 


2 


(i)2. 57982976 


4.1067802 


10" 


2 


1.35219268 


4.9636532 


io- 


2 


(i)2.86633683 


5.0155809 


10" 


2 


1.43415622 


5.2976144 


io- 


2 


(i)3. 15274466 


5.9234204 


10" 


2 


1.51685560 


5.5841641 


io- 


2 


(i)3. 43712011 


6.8307813 


10" 


2 


1.60003382 


5.8244121 


10" 


2 











solution, we use simple polynomial functions. But once 
one numerical solution for a particular equation of state 
is known, Colsys allows the use of these solutions as the 
initial guess. So once a solution is obtained, another so- 
lution can be generated by small modifications of the pa- 
rameters of the previous solution. The smaller the change 
in the parameters, the faster the convergence of the pro- 
grams. The interesting region of the plane (wsr, wj>) can 
be explored by this procedure, looking for null determi- 
nants, in a reasonable amount of time. 
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FIG. 1: Example of Regge- Wheeler function inside the star. 
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FIG. 2: Example of phase function g outside the star. 



RESULTS 



Wc have made several tests on our method using sim- 
pler equations of state. We successfully reproduce data 
from previous works. For example in tableUwe show our 
results for a star of constant density and 2M/R — 0.885, 
in perfect agreement with [46j . 

In Figures Q] and [2] we show an example of calculated 
perturbation functions. These functions correspond to 
a star of 1.4M©, with the high density Sly equation of 
state. The quasi-normal mode is the first axial mode 
wj with frequency v = 8.034:kHz and damping time r = 
29.31/is. In Figure[T]we plot the real and imaginary parts 
of the two interior solutions for Z vs r. In Figure [2] we 
plot the real and imaginary parts of the calculated phase 
function outside the star versus the compactified rotated 
coordinate x. 

In the next sections we will show the results obtained 
for realistic equations of state. 



In this section we present our results for the axial 
quasi-normal modes of neutron stars with realistic equa- 
tions of state. From the equations of state implemented 
using [ij} , we focus the study on the following four equa- 
tions of state: the high density SLy [65(, a equation of 
state for plain nuclear matter; ALF4, a hybrid equa- 
tion of state with mixture of nuclear matter and color- 
flavor locked quark matter (66|; and two equations of 
state for mixed nuclear and hyperon matter, GNH3 67] 
and BGN1H1 [68]. In Figure 0] we show the different 
equations of state in the high density region. For lower 
densities the equation of state used is the Sly [l?} • In 
Figure H] we show the mass radius relation for these equa- 
tions of state. For the study of the quasi-normal modes 
we have considered only stable configurations below the 
maximum mass of each equation of state. The introduc- 
tion of a hyperon core changes the mass radius relation 
and we will see that this also has an imprint on the quasi- 
normal mode spectrum of these configurations. 
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FIG. 3: Density versus pressure in logarithmic scale for the 
four equations of state considered, in the high density region. 
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FIG. 5: Frequency of the fundamental wl mode vs M/R. 
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FIG. 4: Mass vs radius relation for the four equations of state 
considered. 
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FIG. 6: Damping time of the fundamental wl mode vs M/R. 



We start studying the first w modes, in order to ex- 
amine the imprint of each equation of state in the quasi- 
normal modes. 

In Figure [5] we plot the frequency of the fundamental 
wl mode as a function of the compactness of the star, 
M/R, for the different equations of state. We plot from 
compactness 0.1 up to the maximum mass configuration, 
so we consider only stable stars. We observe that for the 
wl modes the presence of hyperon condensates has a clear 
impact on the frequencies of the fundamental w-modes, 
while there is no significant difference between pure nu- 
clear matter and mixed quark- nuclear states. Note how 
the hyperon core is changing the frequency relation for 
high compactness configurations, increasing with M/R 
above a given compactness. In Figure |5] we make a sim- 
ilar plot for the damping time of these configurations. 

In order to use future observations of gravitational 
waves to estimate the mass and the radius of the neutron 



star, as well as to discriminate between different families 
of equations of state, we obtain empirical relations be- 
tween the frequency and damping time of quasi-normal 
modes and the compactness of the star, following (3|.|69j. 
[z3| and [HJ. In Figure [7] we present the frequency of the 
fundamental mode scaled to the radius of each configura- 
tion. It is interesting to note that even with this scaling, 
the softer equations of state that include hyperon matter, 
and in particular, the BGN1H1 equation of state, present 
a quite different behavior than SLy and ALF4 equations 
of state. 

A linear fit can be made to these last two equations of 
state. We fit to 



Lj(Khz) 



1 f M 
R(Km) V R 



D 



(40) 
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TABLE II: Fitting for the wIO modes. Parameters A and B correspond to the linear empirical relation for the frequency (140 \ . 
Parameters a,b and c correspond to the quadratic empirical relation for the damping time (|41|l . 



wIO modes 


SLy 


ALF4 


GNH3 


BGN1H1 


A 


— 148.7 ± 4.5 


-141.4 ± 1.2 


-122.1 ± 1.4 




B 


119.8 ± 1.0 


118.40 ±0.25 


116.67 ±0.25 




x 2 


0.760 


0.039 


0.057 




a 


-1221 ±22 


-1158 ±26 


-1444 ±23 


-1690 ±49 


b 


365.1 ±9.3 


345.3 ±9.8 


425.4 ±8.3 


506 ± 18 


c 


21.63 ±0.89 


23.06 ± 0.88 


17.87 ±0.68 


11.2 ± 1.6 


x 2 


0.050 


0.028 


0.027 


0.084 



TABLE III: Fitting for the wll modes. Parameters A and B correspond to the linear empirical relation for the frequency (|40|l . 
Parameters a,b and c correspond to the quadratic empirical relation for the damping time (|41|l . 



wll modes 


SLy 


ALF4 


GNH3 


BGN1H1 


A 


-689.5 ±5.2 


-688.8 ±6.8 


-736.5 ± 6.8 


-779.4 ± 6.8 


B 


312.5 ± 1.1 


313.5 ± 1.4 


317.2 ± 1.4 


323.3 ± 1.4 


x 2 


0.997 


1.236 


0.372 


1.986 


a 


-1843 ±55 


-1857 ±110 


-2062 ± 25 


-1771 ± 123 


b 


650 ± 23 


654 ± 42 


690.9 ±8.9 


582 ± 45 


c 


14.3 ±2.2 


14.2 ±3.8 


13.10 ±0.73 


21.1 ±3.9 


x 2 


0.311 


0.510 


0.0252 


0.393 



In Table |TT] we present the fitting parameters A,B for 
each one of the equations of state studied. For the Sly 
and for the ALF4 equations very similar results are ob- 
tained. The empirical parameters for SLy and ALF4 are 
compatible with the empirical relation obtained by Ben- 
har et al [7(| for six equations of state. 

For the GNH3 equation of state we still obtain a linear 
relation, but with different parameters as can be seen in 
Table [III The introduction of hyperon matter is affecting 
the empirical relation. 

For the BGN1H1 equation of state the relation of the 
scaled mass with the compactness is not even linear, as 
can be seen in Figure [7] Our results show that the in- 
troduction of hyperon matter changes the dependence of 
the frequency with the compactness. For compactness 
bigger than 0.18, the difference between the equations of 
state with and without hyperon matter is quite impor- 
tant. This feature could be used applying the technique 
from [6t| to test the presence of hyperon matter inside 
neutron stars. Indications in this direction have been 
obtained in a recent work ([HI)- Chatterjee and Bandy- 
opadhyay have studied the wIO mode for other equations 
of state involving hyperons as well as Bose-Einstein con- 
densates of antikaons obtaining that exotic matter can 
be distinguish from nucleons only matter. 

In Figure [5] we present the damping time of the funda- 



mental mode scaled to the mass of each configuration. In 
this case, the results can be fitted to a empirical quadratic 
relation on M/R, as follows: 



1 



1 



-(ps) M(M ) 



(41) 



In Table 11*11 we present the fitting parameters a,b and c. 
Note that for SLy and ALF4 the fits are quite similar, 
and in accordance with the results obtained in ([70]). For 
GNH3 and BGN1H1 equations of state, although the re- 
sult are well fitted to a quadratic relation, the coefficients 
are very different from those obtained for Sly and ALF4 
equations of state. As occurred with the frequency, this 
feature could be used to test the presence of hyperon 
matter in the equation of state. 

Considering the results obtained for the fundamental 
modes of Sly and ALF4 equations of state, it is not pos- 
sible to distinguish between mixed quark-nuclear matter 
composition and pure nuclear matter composition. 

We have developed a similar study with the first ex- 
cited wl-mode. In Figures [jj] and [TU] we plot the scaled 
frequency and damping time for these modes. The same 
empirical relations (|4Q[|41I) . with different parameters can 
be found. In Table [TTTI we present the corresponding em- 
pirical parameters. 



12 





105.0 p 




100.0 - 




95.0 - 


= 






90.0 - 


N 




X 










85.0 - 


3 






80.0 - 




75.0 - 




70.0 - 




a 
3 

3 



0.10 0.13 0.16 0.19 0.22 0.25 0.28 0.31 
M/R 




100.0 



0.10 0.13 0.16 0.19 0.22 
M/R 



0.25 0.28 



0.31 



FIG. 7: Scaled frequency of the fundamental wl mode vs FIG. 9: Scaled frequency of the first excited wl mode vs M/R. 
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FIG. 8: Scaled damping time of the fundamental wl mode vs 
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FIG. 10: Scaled damping time of the first excited wl mode vs 
M/R. 



Concerning the frequency of the wll modes, although 
the linear relation for the scaled frequency is well satisfied 
for all the equations of state, the slope is slightly steeper 
for the GNH3 and BGN1H1. 

Concerning the damping time of the wll modes, the 
situation is quite similar to the frequency. The quadratic 
relation is very well satisfied, but note that the param- 
eters of the quadratic relation are very sensitive to the 
presence of hyperon matter. 

In Figure [TT] we plot the frequency of the fundamental 
wll mode scaled with the mass versus the metric func- 
tion e 2A . It can be seen that for all the equations of state 
studied, there is a minimal compactness, below which the 
wll mode do not exist. The scaled frequency is quite lin- 
ear with e 2A nearby this minimal compactness. Our nu- 
merical scheme allows us to obtain explicitly these limit 
configuration with vanishing real part of the mode with 
a precision of the order of 10 -5 . For the SLy and ALF4 
equations of state, the limit configuration has compact- 



ness M/R = 0.1063. For the GNH3 equation of state the 
limit configuration has compactness M/R = 0.1051, and 
for the BGN1H1, M/R = 0.1061. Note that these con- 
figurations for which the fundamental wll vanishes has a 
compactness quite independent of the equation of state. 
These configurations were conjectured by Wen et al in 
25] . Although in this work they were not able to obtain 
explicity the limit configuration, so the compactness of 
the limit configuration were obtained by extrapolation. 
Our results are in agreement with their extrapolation. 

In Figure [T^] we show the damping time of the wll 
mode scaled with the mass. Although for low compact- 
ness the damping time is quite independent of the equa- 
tion of state and linear in M/R, for compactness above 
0.18 it is very sensitive to the presence of hyperon mat- 
ter, changing completely the linear behavior. The damp- 
ing time does not vanish in the limit compactness as the 
frequency does. By the contrary it tends to a limit r 
dependent of the equation of state. 

We would like to remark the similitude between the 
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SLy and ALF4 frequencies and damping times for every 
neutron star mass. Their first w-modes are practically 
indistinguishable. In fact, it has been shown that for re- 
alistic equations of state like the ALF4, for mixed nuclear 
and quark matter, the mass-radius relation of the stars 
is quite similar to the mass-radius relation of stars com- 
posed of pure nuclear matter 66] . This makes impossible 
to decide if a star presents quark matter in their interior 
by measuring the mass and the radius. Our results show 
that neither the observation of the frequencies or damp- 
ing times of the axial w-mode spectrum of neutron stars 
could allow us to decide if the equation of state presents 
mixed nuclear-quark matter states. On the contrary, the 
determination of the frequencies could tell us if there are 
hyperon condensates inside the star. 

Let us present an empirical relation between the real 
part and the imaginary part of the modes, that seems to 
be independent of the equation of the state. As far as we 
know this has not been obtained before. In Figure [TBI we 
plot, for multiple configurations and all the equations of 



state studied, the real part versus the imaginary part of 
the quasi-normal mode. The modes are in units of the 
central pressure. So to obtain the modes with correct 
units one have to use the following formulas: 



w{Khz) = y/p c (cm 2 ) 



1 



1000 2tt 



r(fis) = (VpcCcm- 2 )— gwi m ) \ 



(42) 



where c is the speed of light in cm/s. It is interesting to 
note the following universal relations: 

For the fundamental wl modes, our results of the imag- 
inary part versus the real part can be fitted to a quadratic 
relation as follows: 

w Jm = (-1.98 ± 0.22)+ (43) 
(0.497 ± 0.030)cjr £ + (0.0193 ± 0.0008)^, 

with x 2 = 0.0624. 

For the first excited wl mode, we obtain the following 
fit: 

w/ m = (-1.33 ±0.22)+ (44) 
(0.495 ± 0.014)wR e + (0.00153 ± 0.00018)^ e , 

with x 2 — 0.118773. It is almost linear with the real 
part. 

For the fundamental wll mode, the real part is 
quadratic with the imaginary part. We obtain the fol- 
lowing fit: 

lu Rc = (4.23 ± 0.27) + (45) 
(0.170 ± 0.024)a; /m + (-0.00524 ± 0.00041)cj| m , 

with x 2 — 0.188. These relations allows to obtain the 
LOi m if U)Re is known, recalling relations of this type in 
other context (Kramers-Kroning relations) . From a prac- 
tical point of view this relation can be used to simplify 
the numerical calculations to obtain quasi-normal modes 
for any equation of state, reducing, in our method, the 
time necessary to find the zeros of the determinant. 

We would like to point out that these empirical re- 
lations can be used to estimate the central pressure of 
the neutron star. For instance, if the frequency and 
the damping time of the fundamental mode of a neu- 
tron star are detected, using (l4*2"j) and (f4*3"|) we are able 
to obtain a value for the central pressure p c . Actually, 
given the precision of the fit, we can only obtain the or- 
der of the magnitude of the central pressure valid for 
any equation of state. If the mass and the radius of the 
neutron star have been measured with another method, 
combining these data with the central pressure, we could 
test the validity of a particular equation of state, using 
the mass versus radius relation for it. In Figure [TJ] we 
show the axial quasi-normal mode spectrum up to the 
fifth excited mode of neutron stars with M = 1AM@ 
for the four equations of state. For stars of this mass 
the spectrum for SLy and ALF4 almost coincide. Hence, 
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these equations of state cannot be distinguished up to the 
fifth mode. The spectrum for equations of state includ- 
ing hyperon matter is slightly different. When the order 
of the mode increases, the frequency and the damping 
time of the modes allow to distinguish between neutron 
stars with and without hyperon matter. So an study of 
different modes could be used to determine which par- 
ticle species are populating the inner layers of the star. 
Finally, we study the effect of core- crust transition pres- 
sure variations in the quasi-normal mode spectrum using 
our model of crust as a surface energy density envelop- 
ing the core. We will consider the SLy equation of state, 
for which the core- crust transition is found at density 
1.29 • I0 14 g/cm 3 . We construct a star of 1.4M with and 
without a surface energy density enveloping the core of 
the star at this core- crust transition region. We can com- 
pare the difference in the fundamental w modes between 
our approximation of crust as a thin shell enveloping the 
core, and the standard calculation. The result is the fol- 
lowing: for the fundamental wl mode, the change in the 



frequency is 0.07% and the variation in the damping time 
is —0.02%. For the fundamental wll mode, the changes 
are bigger. For the frequency 0.3% and for the damp- 
ing time —0.07%. This indicate that this approximation 
makes small changes in the w quasi-normal spectrum. 

In Figures [15] and [16] we present the influence of core- 
crust transition pressure variation in the fundamental w- 
modes. We consider an initial neutron star with IAMq 
and surface energy density at the core- crust transition 
point: 4.8 • 10 32 dyn / cm 2 . Fixing the central pressure, 
we variate the core- crust transition pressure between 
10 32 dyn /cm 2 and 10 33 dyn/cm 2 , and compare with the 
modes for the initial configuration. In Figure [15] we plot 
the relative variation in the frequency. It can be seen 
that the impact of transition pressure variations on the 
frequency is small, between —0.05% and 0.08% for the 
fundamental wl mode, and between —0.24% and 0.35% 
for the fundamental wll mode. In Figure [16] we plot the 
relative variation in the damping time with a similar re- 
sult. Transition pressure variations has a small impact on 
the damping time, between 0.02% and —0.035% for the 
fundamental wl mode, and between 0.06% and —0.1% 
for the fundamental wll mode. 



VI. SUMMARY 

In this paper we have studied the axial w quasi-normal 
modes for several equations of state. In particular, we 
consider the influence of the presence of hyperons in the 
core of the neutron star. We have obtained that, in prin- 
ciple, for enough compact stars we can distinguish be- 
tween stars with hyperons and stars without hyperons in 
their core. 

In order to do that we have developed a new algo- 
rithm, and a new numerical code to calculate w-axial 
quasi-normal modes for neutron stars with realistic equa- 
tions of state. 

In the exterior of the star we use the phase of the 
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Regge-Whcclcr function, which is not an oscillatory func- 
tion as the Regge- Wheeler function. In order to generate 
properly outgoing quasi-normal waves, removing the di- 
vergence towards infinity, we have used Exterior Complex 
Scaling, allowing us to impose constringent enough con- 
ditions at infinity by the compactification of the exterior 
region. This treatment allows us to grant that the exte- 
rior solution corresponds to an outgoing wave, without 
any possibility of incoming wave contamination. What 
is new in this part of the algorithm is that we allow for 
a variable angle in the complex scaling, which makes it 
possible to obtain new modes and better precision than 
in previous works. 

In the interior of the star we obtain the general solution 
of the Regge- Wheeler equation, as a linear combination 
of two numerically generated independent solutions for 
neutron star with realistic equations of state and reg- 
ularity conditions in the center of the star. Then the 
interior and the exterior solution have to be matched ad- 
equately in the border of the star. Only for some values 
of the frequency ui and damping time r it is possible 
to satisfy these multi-boundary conditions: conditions in 
the center of the star, in the border of the star and at 
the spatial infinity. Our algorithm reduces the search for 
these values of ui and r to the search of the zeros of the 
determinant of a matrix evaluated in the border of the 
star. 

A complete study of the junction conditions in a sur- 
face of constant pressure (also including the case p = 0, 
the border of the star) has been done using intrinsic Dar- 
mois formulation. We include the general case with a 
surface layer in the boundary that can be used to study 
situations with a rapid variation in the density (for in- 
stance, in the crust-core transition region). 

It is interesting to note some technical points of the 
algorithm: we have used Colsys, a collocation package, 
to integrate the differential equations given its ability to 
treat multi-boundary conditions and it has great flexi- 



bility available to impose the precision on the different 
functions. We integrate the background (zero-order) and 
the perturbation (first-order) at the same time to ob- 
tain the precision in the zero order needed to calculate 
the perturbation with the precision required. Usually for 
high excited modes the Regge- Wheeler function presents 
an important number of oscillations inside the star. More 
precision on the static functions (zero order) is needed, 
but as the functions are all resolved at once, the pro- 
grams just adapt the mesh until the convergence criteria 
is achieved. In order to be able to integrate the zero order 
and the perturbation at the same time, we have added a 
trivial equation for the border of the star R. Also, it is 
interesting to note, that the continuation method is used 
to generate families of solutions varying slowly a given 
parameter, and that a large region of the complex plane 
of lo has been explored looking for zeros of the determi- 
nant condition. 

One expects that the spectrum of quasi-normal modes 
has to depend on the equation of the state of the neutron 
star. Hence, the observation of the spectrum of quasi- 
normal modes will impose constrains on the equation of 
the state (EOS). To do that, in an effective manner, it 
is interesting to have a parameterization of the EOS in 
terms of a small number of parameters. In this work 
we have used a piece-wise polytropic approximation by 
Jocelyn et al. |17| . which fits, with a very good precision, 
a large family of EOS. 

To study the influence of the presence of hyperons in 
the core of the neutron star we concentrate our study in 
four EOS: SLy (nuclear matter), ALF4 (mixed nuclear 
and quark matter), GNH3 and BGN1H1 (with hyperon 
condensates) . 

We have obtained phenomenological relation for the 
frequency and damping time with M/R for wl and wll 
modes indicating that we can not differentiate the results 
for SLy(nuclear matter) and ALF4 (mixed nuclear and 
quark matter). But the presence of hyperons in GNH3 
and BGN1H1 EOS influence enough the spectrum to 
make possible to distinguish the spectrum of a star with 
hyperons and a star without hyperons if the compactness 
is large enough. 

Other interesting result is that we have obtained a new 
phenomenological relation between the real part and the 
imaginary part of the frequency of the w quasi-normal 
modes valid for the four EOS. This relation allows to 
obtain an estimation of the central pressure of the star 
independent of the EOS, if the frequency and the damp- 
ing time of the w fundamental mode of a neutron star is 
detected. 

Also, the precision of our algorithm allows us to con- 
struct explicitly the universal low compactness limiting 
configuration around M/R = 0.106 for which the funda- 
mental wll mode vanishes. The existence of this config- 
uration had been conjectured by (25|. 

Finally, we have studied the influence of changes in 
the core- crust transition pressure obtaining that the in- 
fluence is very small. The relative changes in the fre- 
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quency for a variation of the core- crust transition from 
I0 32 dyn/cm 2 to 10 33 dyn/cm 2 are between 0.02% and 
—0.035% for the fundamental wl mode, and between 
0.06% and —0.1% for the fundamental wll mode. In 
the damping time the relative changes are between 0.02% 
and —0.035% for the fundamental wl mode, and between 
0.06% and —0.1% for the fundamental wll mode. 



tron stars. J.L.B wants to thank also his kind hospi- 
tality. Also, we are very grateful to Prof. Gabriel A. 
Galindo for his help concerning the Exterior Complex 
Scaling method. This work was supported by the Span- 
ish Ministerio de Ciencia e Innovacion, research project 
FIS2011-28013. J.L.B was supported by the Spanish Uni- 
versidad Complutense de Madrid. 



Acknowledgments 

We thank Prof. Kostas D. Kokkotas for very help- 
ful discussions concerning quasi-normal modes of neu- 



[20] 



M. Pitkin, S. Reid, S. Rowan, and J. Hough, Living Re- 
views in Relativity 14 (2011). 

B. Sathyaprakash and B. F. Schutz, Living Reviews in 
Relativity 12 (2009). 

N. Andersson, V. Ferrari, D. Jones, K. Kokkotas, B. Kr- 
ishnan, J. Read, L. Rezzolla, and B. Zink, General Rela- 
tivity and Gravitation 43, 409 (2011). 
N. Andersson and K. D. Kokkotas, Physical Review Let- 
ters 77, 4134 (1996). 

M. Shibata, K. Taniguchi, and K. b. o. Uryfi, Phys. Rev. 
D 71, 084021 (2005). 

M. Shibata and K. Uryu, Progress of Theoretical Physics 
107, 265 (2002). 

R. Oechslin, H.-T. Janka, and A. Marek, A&A 467, 395 
(2007). 

T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, 
Phys. Rev. D 81, 123016 (2010). 

J. S. Read, C. Markakis, M. Shibata, K. b. o. Uryu, 
J. D. E. Creighton, and J. L. Friedman, Phys. Rev. D 
79, 124033 (2009). 

J. Abadie, B. P. Abbott, R. Abbott, R. Adhikari, 
P. Ajith, B. Allen, G. Allen, E. Amador Ceron, R. S. 
Amin, S. B. Anderson, et al. (The LIGO Scientific Col- 
laboration), Phys. Rev. D 83, 042001 (2011). 
K. D. Kokkotas and B. Schmidt, Living Reviews in Rel- 
ativity 2 (1999). 

H.-P. Nollert, Classical and Quantum Gravity 16, R159 
(1999). 

L. Rezzolla, Black Holes and Relativistic Stars p. 47 
(2003). 

A. Haensel, P. Potekhin and D. Yakovlev, Neutron stars: 
Equation of state and structure, Astrophysics and space 
science library (Springer, 2007). 

N. Glendenning, Compact stars: nuclear physics, parti- 
cle physics, and general relativity, Astronomy and astro- 
physics library (Springer, 2000). 

H. Heiselberg and M. Hjorth- Jensen, Physics Reports 
328, 237 (2000). 

J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, 

Phys. Rev. D 79, 124032 (2009). 

L. Lindblom, Phys. Rev. D 82, 103011 (2010). 

O. Benhar, E. Berti, and V. Ferrari, Monthly Notices of 

the Royal Astronomical Society 310, 9 (1999). 

K. D. Kokkotas, T. A. Apostolatos, and N. Andersson, 

Monthly Notices of the Royal Astronomical Society 320, 



307 (2001). 

[21] O. Benhar, V. Ferrari, and L. Gualtieri, Physical Review 

D 70, 124015 (2004). 
[22] O. Benhar, Mod.Phys.Lett. A20, 2335 (2005). 
[23] V. Ferrari and L. Gualtieri, Gen.Rel.Grav. 40, 945 

(2008). 

[24] D. Chatterjee and D. Bandyopadhyay, Phys. Rev. D80, 
023011 (2009). 

[25] D.-H. Wen, B.-A. Li, and P. G. Krastev, Phys.Rev. C80, 
025801 (2009). 

[26] T. Regge and J. A. Wheeler, Physical Review 108, 1063 
(1957). 

[27] F. J. Zerilli, Phys. Rev. Lett. 24, 737 (1970). 
[28] K. Thorne and A. Campolattaro, Astrophys. J. 149, 591 
(1967). 

[29] R. Price and K. Thorne, Astrophys. J. 155, 163 (1969). 
[30] K. Thorne, Astrophys. J. 158, 1 (1969). 
[31] K. Thorne, Astrophys. J. 158, 997 (1969). 
[32] A. Campolattaro and K. Thorne, Astrophys. J. 159, 847 
(1970). 

[33] L. Lindblom and S. Detweiler, Astrophys. J. Supplement 

Series 53, 73 (1983). 
[34] S. Detweiler and L. Lindblom, Astrophys. J. 292, 12 

(1985). 

[35] S. Chandrasekhar and V. Ferrari, Proceedings of the 
Royal Society of London. Series A: Mathematical and 
Physical Sciences 432, 247 (1991). 

[36] S. Chandrasekhar, V. Ferrari, and R. Winston, Proceed- 
ings of the Royal Society of London. Series A: Mathe- 
matical and Physical Sciences 434, 635 (1991). 

[37] S. Chandrasekhar and V. Ferrari, Proceedings of the 
Royal Society of London. Series A: Mathematical and 
Physical Sciences 434, 449 (1991). 

[38] Y. Kojima, Phys. Rev. D 46, 4289 (1992). 

[39] Y. Kojima, N. Andersson, and K. D. Kokkotas (1995). 

[40] K. D. Kokkotas, Monthly Notices of the Royal Astro- 
nomical Society 268, 1015 (1994). 

[41] K. D. Kokkotas and B. F. Schutz, Monthly Notices of the 
Royal Astronomical Society 255, 119 (1992). 

[42] N. Andersson, K. D. Kokkotas, and B. F. Schutz, 
Monthly Notices of the Royal Astronomical Society 274, 
9 (1995). 

[43] S. Chandrasekhar and S. Detweiler, Proceedings of the 
Royal Society of London. A. Mathematical and Physical 
Sciences 344, 441 (1975). 



17 



[44] N. Andersson, Y. Kojima, and K. D. Kokkotas, Astro- 
physical Journal 462, 855 (1996). 

[45] N. Andersson, Proceedings of the Royal Society of Lon- 
don. Series A: Mathematical and Physical Sciences 439, 
47 (1992). 

[46] L. Samuelsson, N. Andersson, and A. Maniopoulou, Clas- 
sical and Quantum Gravity 24, 4147 (2007). 

[47] E. W. Leaver, Proceedings of the Royal Society of Lon- 
don. Series A: Mathematical and Physical Sciences 402, 
285 (1985). 

[48] M. Leins, H. P. Nollert, and M. H. Soffel, Phys. Rev. D 

48, 3467 (1993). 
[49] G. Darmois, Memorial des Sciences Mathematiqucs 

(Paris:Gauthier-Villars) 25 (1927). 
[50] A. Lichnerowicz, Theories Relativistes de la Gravitation 

et de UElectromagnetism (Maison, Paris, 1955). 
[51] S. O'Brien and J. Syngr, Comm. Dublin Inst. Adv. Stud. 

Ser. A 9, 1 (1952). 
[52] W. Israel, Nuovo Cimento B B44, 1 (1966). 
[53] M. Mars and J. M. Senovilla, Class.Quant.Grav. 10, 1865 

(1993). 

[54] R. Vera, Class.Quant.Grav. 19, 5249 (2002). 

[55] M. Mars, Class.Quant.Grav. 22, 3325 (2005). 

[56] L. M. Gonzalez-Romero, Phys. Rev. D 67, 064011 (2003). 

[57] M. A. H. MacCallum, M. Mars, and R. Vera, Phys. Rev. 



D 75, 024017 (2007). 

[58] H. Sotani, K. Tominaga, and K.-i. Maeda, Phys. Rev. D 
65, 024010 (2001). 

[59] C. Misner, K. Thorne, and J. Wheeler, Gravitation, 
Physics Series (W. H. Freeman, 1973). 

[60] U. Ascher, J. Christiansen, and R. D. Russell, Mathe- 
matics of Computation 33, pp. 659 (1979). 

[61] J. Aguilar and J. Combes, Commun. Math. Phys. 22, 
269 (1971). 

[62] E. Balslev and J. Combes, Commun. Math. Phys. 22, 
280 (1971). 

[63] B. Simon, Commun. Math. Phys. 27, 1 (1972). 

[64] G. Alvarez, R. Damburg, and H. Silverstone, Phys. Rev. 

A 44, 3060 (1991). 
[65] F. Douchin and P. Haensel, A&A 380, 151 (2001). 
[66] M. Alford, M. Braby, M. Paris, and S. Reddy, Astrophys. 

J. 629, 969 (2005). 
[67] N. Glendenning, Astrophys. J. 293, 470 (1985). 
[68] S. Balberg and A. Gal, Nuclear Physics A 625, 435 

(1997). 

[69] N. Andersson and K. D. Kokkotas, Monthly Notices of 
the Royal Astronomical Society 299, 1059 (1998). 

[70] O. Benhar, E. Berti, and V. Ferrari, Monthly Notices of 
the Royal Astronomical Society 310, 797 (1999). 



